
\chapter{\label{sec:Online-Optimization-of}Real-time Optimization of DAC }

Distributed average consensus (DAC) algorithm is widely used in many
applications. It utilizes matrix iteration to find the dominant eigenvector.
To minimize the required number of iterations, the algorithm needs
to be optimized. However, this optimization needs the knowledge of
network topology, which is very hard to obtain for an individual agent
in distributed networks.  Thus, optimal step length and forgetting
factor need to be calculated offline and forwarded to every agent.
 To solve this problem, we proposed a distributed real-time optimization
technique so that each node can estimate these optimal parameters
individually. In addition, the method is based on constant first-order
DAC itself, so it will not stop the consensus process. The result
shows that a numerical error due to quantization would exist in the
distributed solution. It will increase as the network becomes larger.
Thus, a numerical  technique is introduced  to mitigate the error.
The estimated parameters after mitigation do not obviously  decline
the performance of higher-order DAC when network size is smaller than
a threshold. 




\section{Introduction}

In many applications, such as time synchronization \cite{Schenato2011},
cooperative control of vehicles \cite{Yang2010}, formation control
\cite{Olfati-Saber2012} and WSNs \cite{Hlinka2012}, it is often
necessary that a group of agents in a distributed system can agree
on certain quantities. An example application is distributed detection
of a moving target by wireless sensor networks. Suppose each sensor
is observing the target coordinates but the output is corrupted by
independent and identically distributed zero-mean Gaussian noise,
to minimize the interference from the noise, the sensors need to take
the average of all initial values. 

The problem of how to achieve this average in a distributed system
is called the \textit{average consensus problem},  which is solved
by distributed average consensus (DAC) algorithms. 

When the moving target is highly dynamic or the sensors need to sample
at a very high frequency, it requires that  the DAC algorithm returns
the result in a short time. Thus, many efforts have been devoted to
optimize the algorithm.

The DAC algorithms can be divided into asymptotic and non-asymptotic
algorithms. Asymptotic algorithms have been proved to be robust against
 topology changes and they play important roles in practice \cite{Ren2007}.
 The optimization of asymptotic algorithms is to minimize the sub-dominate
eigenvalue of a weight matrix \cite{Asensio-Marco2012}\cite{Xiao2004}\cite{Xiong2010}.
 In addition, \cite{Kokiopoulou2007} deals with the iteration acceleration.
However, these optimization of DAC algorithms are centralized methods,
which means a centralized node calculates the optimal parameters and
forwards them to the whole network \cite{Xiong2010}. 

To enable the whole system work distributively, the optimization should
be distributed and the DAC algorithms should be robust against topology
changes. A distributed method inspired by the gossip algorithm \cite{Boyd2006}
can be used to optimize the first order DAC but it converges very
slowly.   The method involves triple nested distributed matrix iterations.
The inner iteration has to converge to a certain range so that the
iteration outside can return the right result. Thus, It is not surprising
that it could not finish in a reasonable time when the network size
is large.

For non-asymptotic DAC algorithms, such as finite-time \cite{Sundaram2007}
and adaptive filter DAC algorithm \cite{Cavalcante2010}, they find
the average in some sophisticated ways. Sumdaram and Hadjicostis \cite{Sundaram2007}
verify that there exists a filter that can estimate the consensus
value. Cavalcante and Mulgrew \cite{Cavalcante2010} follow Sundaram
and Hadjicostis's work to propose an adaptive algorithm to find the
filter. The optimization for both of them is to minimize the number
of necessary iteration before a FIR filter is estimated. However,
they are not robust against topology changes. Both of them have a
first-order DAC running in the background and the local values over
time are taken as inputs of the filter estimation algorithm. As a
result, if the network topology changes, these algorithms have to
be terminated, as outdated information during the filter estimation
will lead to a wrong answer.



After investigating these problems, we intend to find a distributed
optimization method for the constant first-order DAC and higher-order
DAC algorithms. Because in centralized optimization methods, optimal
parameters of these algorithms are only related to the eigenvalues
of Laplacian matrix of the network \cite{Xiong2010}, if we could
estimate these eigenvalues in a distributed manner, these centralized
methods could be carried out distributively. 

Consequently, a distributed eigenvalue estimation algorithm is proposed
in this chapter. In contrast to other distributed algorithms \cite{Kempe2008}\cite{Franceschelli2009}\cite{Yang2010},
initialization of the proposed algorithm is actually the first-order
DAC itself. Therefore, first-order DAC will not be interrupted during
the optimization and algorithm complexity and communication cost can
be dramatically reduced.  

However, the distributed solution has a numerical error due to quantization,
which may decline the algorithm performance. Therefore, a least mean
square solution is obtained to mitigate the numerical error. When
using the floating point number in double format and the network size
is smaller than 32, the numerical error after mitigation does not
dramatically decline the performance and the proposed method is applicable. 





The rest of this chapter is structured as follows. First,  the distributed
real-time optimization of DAC will be given in section \ref{sec:distributed-Eigenvalue-Estimati}.
Second, the mitigation of numerical error  will be proposed in section
\ref{sec:Mitigation-of-Numerical}. Third, the algorithm complexity
will be analysed in section \ref{sec:Algorithm-Complexity}. Fourth,
in section \ref{sec:Simulation-of-Eigenvalue}, the performance of
DAC using the distributed real-time optimization  will be  analysed
and compared with the centralized one. Finally, the conclusion will
be given in section \ref{sub:Conclusion}.


\section{\label{sec:distributed-Eigenvalue-Estimati}Real-time optimization
of DAC }

Traditional optimization of HO-DAC or CFO-DAC requires a centralized
node to gather information of the Laplacian matrix \cite{Xiong2010}.
Without the spectrum of Laplacian matrix, each node could only choose
a non-optimal point $\left(\epsilon,\gamma\right)$ in the boundary
of  the convergence region.

Because in section \ref{sub:Discrete-High-Order}, it is shown that
the optimal parameters $\epsilon_{opt},\gamma_{opt}$ of HO-DAC are
only related to $\lambda_{i}\left(L\right)$, to enable the distributed
optimization, the key is to estimate these eigenvalues in a distributed
manner.  

In fact, there are some decentralized techniques \cite{Kempe2008}\cite{Franceschelli2009}\cite{Yang2010}
to estimate the eigenvalues. However, they are not designed for DAC
algorithm and will involve complex and costly initialization. In addition,
as they are also matrix iteration algorithms, DAC algorithm has to
be stopped during the eigenvalue estimation. 

In contrast, the distributed real-time optimization and CFO-DAC can
be running simultaneously  and algorithm complexity and communication
cost can be dramatically reduced. Because the initialization of the
proposed algorithm is to store a number of local values  obtained
by  CFO-DAC, the distributed system can still work using non-optimal
parameters at the very beginning just after the deployment. After
a number of iterations of CFO-DAC, these eigenvalues could be estimated
and better parameters could be used in the next iterations.

  

Because $\lambda_{i}\left(W\right)$ is the root of the characteristic
polynomial $p(\lambda)=\prod_{i=1}^{m}\left(\lambda-\lambda_{i}\right)^{r_{i}}=\lambda^{D}+a_{D-1}\lambda^{D-1}+\ldots+a_{0}=0,$
the distributed eigenvalues estimation can be cast into calculating
the coefficients $\left\{ a_{j}\right\} $. 

To avoid complicated analysis, we assume duplex communication is possible
in each link. Therefore, the  weight matrix $W$ is symmetric, diagonalizable
and has $n$ linear independent eigenvectors.  To calculate $\left\{ a_{j}\right\} $,
we need to use the following theorem:
\begin{thm}
\label{thm:LV predictor a_i-1} Suppose an undirected graph ${\cal G}$
with associated weight matrix $W\in\mathbb{R}^{n\times n}$ which
satisfies the conditions in theorem \ref{thm:convergence condition}.
For the DAC iteration $\mathbf{x}(k+1)=W\mathbf{x}(k)$, local value
\textup{$x_{i}\left(k\right)$}  is equal to a linear combination
of  itself in $D$ previous time steps, i.e. $x_{i}\left(k\right)=-a_{D-1}x_{i}\left(k-1\right)-\ldots-a_{1}x_{i}\left(k-D+1\right)-a_{0}x_{i}\left(k-D\right)$,
where $k\geq D$ and $D$ is a certain number.\end{thm}
\begin{proof}
Since all the eigenvectors of $W$ consist a basis of $\mathbb{R}^{n}$,
any initial value vector can be decomposed into a linear combination
of these eigenvectors, $\mathbf{x}\left(0\right)=\alpha_{1}\mathbf{e}_{1}+\alpha_{2}\mathbf{e}_{2}+\ldots+\alpha_{n}\mathbf{e}_{n}$.
For $k=1,2,3,...$ we have
\begin{eqnarray}
\mathbf{x}\left(k\right) & = & W^{k}\mathbf{x}\left(0\right)=\sum_{i=1}^{n}\alpha_{i}\lambda_{i}^{k}\mathbf{e}_{i}\label{eq:x(k) decomposition-1}
\end{eqnarray}
where $\lambda_{i}$ is the eigenvalue of $W$.  Substitute   \prettyref{eq:x(k) decomposition-1}
into $p\left(\lambda\right)$, we have
\begin{equation}
\mathbf{x}\left(k\right)+a_{D-1}\mathbf{x}\left(k-1\right)+\ldots+a_{0}\mathbf{x}\left(k-D\right)=\mathbf{0}_{n\times1}\label{eq:x(k+D) vector predictor}
\end{equation}
 $\forall v_{i}\in{\cal V}$, its local value $x_{i}\left(k\right)$
is the $i^{th}$ component of $\mathbf{x}\left(k\right)$. After a
little evolution, we have
\begin{equation}
x_{i}\left(k\right)=-a_{D-1}x_{i}\left(k-1\right)-\ldots-a_{0}x_{i}\left(k-D\right).\label{eq:x(k+D) predictor}
\end{equation}
In other words, $x_{i}\left(k\right)$ can be predicted by a finite
impulse response filter, if $k\geq D$. \end{proof}
\begin{lem}
The sum of coefficients \textup{$a_{D-1},\ldots,a_{1},a_{0}$} is
equal to $-1$\textup{. }\end{lem}
\begin{proof}
 Since  \prettyref{eq:x(k+D) vector predictor} satisfied as $k\to\infty$,
we have $\bar{\mathbf{x}}+a_{D-1}\bar{\mathbf{x}}+\ldots+a_{1}\bar{\mathbf{x}}+a_{0}\bar{\mathbf{x}}=\mathbf{0}$.
Then, cancelling the vector $\bar{\mathbf{x}}$ will obtain $\sum_{j=0}^{D-1}a_{j}=-1$. 
\end{proof}

\subsection{Find the Polynomial Coefficients}

If more local values are available, node $v_{i}$ could list a number
of equations similar to \prettyref{eq:x(k+D) predictor}. Once they
are sufficient to construct a matrix, we can have an important result.

Define the function $\mathbf{y}_{i}=\mathbf{y}_{i}\left(k,D_{i}\right)=\left[x_{i}\left(k+1\right),x_{i}\left(k+2\right),\ldots,x_{i}\left(k+D_{i}\right)\right]^{T}\in\mathbb{R}^{D_{i}}$
which outputs a \textit{local value history vector} of $v_{i}$ from
time $k+1$ to $k+D_{i}$. Then, define a function $T_{i}=T_{i}\left(k,D_{i}\right)\in\mathbb{R}^{D_{i}\times D_{i}}$
that outputs a Toeplitz matrix with $x_{i}\left(k\right)$ on the
diagonal.  $T_{i}\left(k,D_{i}\right)=$
\begin{align}
\left[\begin{array}{cccc}
x_{i}\left(k\right) & x_{i}\left(k-1\right) & \ldots & x_{i}\left(k-D_{i}+1\right)\\
x_{i}\left(k+1\right) & x_{i}\left(k\right) & \cdots & x_{i}\left(k-D_{i}+2\right)\\
\vdots & \vdots & \ddots & \vdots\\
x_{i}\left(k+D_{i}-1\right) & x_{i}\left(k+D_{i}-2\right) & \cdots & x_{i}\left(k\right)
\end{array}\right]\label{eq:Toeplitz_i simple}
\end{align}
Besides, let $\mathbf{a}_{i}\left(D_{i}\right)=\left[a_{i,D_{i}-1},\ldots,a_{i,1},a_{i,0}\right]^{T}\in\mathbb{R}^{D_{i}}$
be a vector  to store the coefficients calculated at $v_{i}$. Finally,
we have the solution of $\mathbf{a}_{i}\left(D_{i}\right)$ given
by
\begin{equation}
\mathbf{a}_{i}\left(D_{i}\right)=-T_{i}^{-1}\left(k,D_{i}\right)\mathbf{y}_{i}\left(k,D_{i}\right).\label{eq:Toeplitz Eq.}
\end{equation}


Toeplitz matrix is a special type of matrix and can be inverted by
Levinson's algorithms in the polynomial time of order $D_{i}^{2}$,
rather than the order of $D_{i}^{3}$ in general case (for example
by LU decomposition) \cite{Prass2007}.     


\subsection{Estimated Eigenvalues of Laplacian Matrix}

After node $v_{i}$ could calculate the coefficients vector $\mathbf{a}_{i}\left(D_{i}\right)$,
it will construct a local polynomial $p_{i}\left(\lambda\right)$
and find the roots of the polynomial. Then, the local eigenvalues
spectrum of $W$ at $v_{i}$ is obtained, which is defined by $\hat{S}_{i}\left(W\right)=\left\{ \hat{\lambda}_{j}^{\left(i\right)}\left(W\right)\right\} =\left\{ \lambda|p_{i}\left(\lambda\right)=0\right\} $,
 $j=1,\ldots,D_{i}$. Because $W=I_{n}-\epsilon L$, an estimated
Laplacian spectrum at $v_{i}$ could be obtained, which is denoted
by $\hat{S}_{i}\left(L\right)=\left\{ \hat{\lambda}_{j}^{\left(i\right)}\left(L\right)\right\} $,
where $\hat{\lambda}_{j}^{\left(i\right)}\left(L\right)=\frac{1-\hat{\lambda}_{j}^{\left(i\right)}\left(W\right)}{\epsilon},\ j=1,2.\ldots,D_{i}$.


\subsection{Eigenvalues missing in local spectrum }

In simulation, some eigenvalues are missed in some of the local eigenvalues
spectrums. The reason is because sometimes the Toeplitz matrix $T_{i}$
with the original size $D$ will loss rank, and node $v_{i}$ could
only build a smaller Toeplitz matrix with size $D_{i}$ so that the
solution is unique. As a result, the local polynomial $p_{i}\left(\lambda\right)$
will reduce to $D_{i}$ ($D_{i}\leq D\leq n$) coefficients. In addition,
$D_{i}$ may different from one to another.

The example of a node $v_{i}$ who can not estimate a eigenvalue $\lambda_{l}$
is as follows. Suppose $\mathbf{x}\left(k\right)=W^{k}\mathbf{x}\left(0\right)=\sum_{j=1}^{n}\alpha_{j}\lambda_{j}^{k}\mathbf{e}_{j}$,
if for any eigenvalue $\lambda_{l}=\lambda_{l}\left(W\right)$, the
associated eigenvector has $i'th$ component equals to zero, i.e $\mathbf{u}_{d}^{T}\mathbf{e}_{l}=0$,
where $\mathbf{u}_{i}$ is a all-zero vector except the $i'th$ component
is one, then the local value at the node $v_{i}$, denoted by $x_{i}\left(k\right)=\mathbf{u}_{i}^{T}\sum_{j=1}^{n}\alpha_{j}\lambda_{j}^{k}\mathbf{e}_{j}$,
will not contain any information of $\lambda_{l}$. 

However, the following theorem in \cite{Sundaram2007} assure that
the roots of local polynomial is the same roots of minimal polynomial
of $W$. 
\begin{thm}
The local polynomial $p_{i}\left(\lambda\right)$ of node $v_{i}$
divides the minimal polynomial $p\left(\lambda\right)$ of $W$ for
all $1\leq i\leq n$. (for proof, see \cite{Sundaram2007})
\end{thm}
For this reason, every node need to exchange its estimated eigenvalue
spectrum with its neighbors in order to recover the full eigenvalue
spectrum of $W$. To assure that no eigenvalue is missed in the final
eigenvalue spectrum, for any eigenvalue $\lambda_{l}$$\left(W\right)$,
it must be estimated by at least one node in the network. 
\begin{thm}
Suppose a graph ${\cal G}$ whose associated weight matrix $W$ satisfies
the convergence condition, and initial value vector $\mathbf{x}\left(0\right)\in\mathbb{R}^{n}$
is chosen randomly. If all nodes estimate the eigenvalue of $W$ by
the proposed method using local values that are available, any eigenvalue
$\lambda_{l}\left(W\right)$ could be estimated by at least one node
in the network.\end{thm}
\begin{proof}
For the case that matrix $W$ is diagnosable, the proof can be easier.
Since $W^{k}=P\Lambda^{k}P^{-1}=\sum_{i=1}^{n}\lambda_{i}^{k}\mathbf{e}_{i}\mathbf{e}_{i}^{T}$,
we have $\mathbf{x}\left(k\right)=W^{k}\mathbf{x}\left(0\right)=\sum_{i=1}^{n}\alpha_{i}\lambda_{i}^{k}\mathbf{e}_{i}$,
and $\alpha_{i}\neq0$, for all $i$. If an node $v_{d}$ could not
estimation a eigenvalue $\lambda_{l}=\lambda_{l}\left(W\right)$,
which means the information of $\lambda_{l}$ is missed in $x_{d}(k)$.
This happens if and only if the $d^{th}$ component of $\mathbf{e}_{l}$
is zero, or the linear combination of eigenvectors associated with
$\lambda_{l}$ is zero component at $d$. If no node in the network
could estimate $\lambda_{l}$, the associated eigenvector or the linear
combination of associated eigenvectors is $\mathbf{0}_{n\times1}$.
However, it is not possible for a diagnosable matrix. Since diagnosable
matrix has $n$ independent eigenvectors, the linear combination of
eigenvectors can not be equal to $\mathbf{0}_{n\times1}$ as well.
Therefore, there exist at least one node in the network could estimate
$\lambda_{l}$.

The proof can be generalized to non-diagnosable matrix with the help
of Jordan decomposition of $W$. 
\begin{eqnarray*}
W^{k} & = & UJ^{k}U^{-1}\\
 & = & U\left[\begin{array}{ccc}
J_{1}^{k}\\
 & \ddots\\
 &  & J_{\hat{m}}^{k}
\end{array}\right]U^{-1}
\end{eqnarray*}
where $\hat{m}$ is the number of Jordan blocks. Therefore, $\mathbf{x}\left(k\right)=UJ^{k}U^{-1}\mathbf{x}\left(0\right)$.
\end{proof}
Assume that all node miss an eigenvalue $\lambda_{l}$ in their estimated
eigenvalue spectrums, $l=1,2,\ldots,\hat{m}$. This means all the
terms contain $\lambda_{l}$ are multiplied by zero. Since $\mathbf{x}\left(0\right)$
is chosen randomly in $\mathbb{R}^{n}$, we can only have $W^{k}=UJ^{k}U^{-1}$
with all terms involving $\lambda_{l}$ are equal to zero, thus 
\begin{equation}
U\left[\begin{array}{ccc}
\mathbf{0}\\
 & J_{l}^{k}\\
 &  & \mathbf{0}
\end{array}\right]U^{-1}=\mathbf{0}_{n\times n}\label{eq:Jordan decom.under assumption}
\end{equation}
Since $J_{l}^{k}\neq\mathbf{0}$ and the matrix $U$ is consist by
linear independent vectors, the above equation \prettyref{eq:Jordan decom.under assumption}
can not be valid. Thus, the assumption is not true and $\lambda_{l}\left(W\right)$
could be estimated by at least one nodes in the network.


\section{\label{sec:Mitigation-of-Numerical}Mitigation of Numerical Error
of Eigenvalues }

Due to the limited accuracy of the floating point number, the solution
obtained by the proposed method has a numerical error. To mitigate
the  error, one of the options is to use floating number with more
bits. However, communication cost is increased in each iteration. 

As floating point number in double  format is  available on most computers,
we apply our numerical mitigation  on this basis.

The idea is to have more equations similar to  \prettyref{eq:x(k+D) predictor}
and involve the Moore-Penrose pseudo-inverse to find the least mean
square solution. Thus, Toeplitz matrix in  \prettyref{eq:Toeplitz Eq.}
should be replaced by a matrix constructed by some other way, whose
height is larger than width. The local value history vector $\mathbf{y}_{i}\left(k,D\right)$
on the left of  \prettyref{eq:Toeplitz Eq.} is also expanded accordingly. 

There are two ways to expand the matrix. First, the new matrix can
be built by concatenating   Toeplitz matrix $T_{i}$  and $T_{j},v_{j}\in{\cal N}_{i}$
along the column.  As $x_{j}$ is available for node $v_{i}$, this
improvement will not increase the communication cost. 

Second, more than one instances of CFO-DAC algorithms could be carried
out to obtain more useful local samples. Let $\mathbf{x}_{1}\left(0\right),\mathbf{x}_{2}\left(0\right),\ldots,\mathbf{x}_{N}\left(0\right)$
denote $N$ different and independent initial local value vectors.
 Each one of them is used to reinitialize an instance of CFO-DAC and
  each instance of CFO-DAC is iterated for at least $2n+2$ steps.
During this initialization, each node $v_{i}$ will store the local
values obtained by each instance of CFO-DAC.

To construct the new matrix, first let $T_{i,s}=T_{i,s}\left(D_{i}-1,D_{i}\right)$
be the Toeplitz matrix of $v_{i}$ at $s$ instance of CFO-DAC, with
$x_{i,s}\left(D_{i}-1\right)$ on the diagonal, where $x_{i,s}$ is
the local value of node $v_{i}$ at $s$ instance of CFO-DAC. Second,
concatenating $T_{i,s}$ and $T_{j,s},\forall v_{j}\in{\cal N}_{i}$
will obtain $M_{i,s}$.
\begin{equation}
M_{i,s}=\left[T_{i,s}^{T},T_{j_{1},s}^{T},\ldots,T_{j_{\left|{\cal N}_{i}\right|},s}^{T}\right]^{T}
\end{equation}
where $s=1,...,N$.  Finally, concatenating the matrix $M_{i,s}$
will construct another matrix, 
\begin{equation}
\tilde{M}_{i}=\left[M_{i,1}^{T},M_{i,2}^{T},\ldots,M_{i,N}^{T}\right]^{T}.
\end{equation}
The width of $\tilde{M}_{i,s}$, denoted by $D_{i}$, is the maximum
integer so that the matrix $\tilde{M}_{i,s}$ has full column rank. 

On the other hand, let $\mathbf{y}_{i,s}=\mathbf{y}_{i,s}\left(D_{i}-1,D_{i}\right)=\left[x_{i,s}\left(D_{i}\right),x_{i,s}\left(D_{i}+1\right),\dots,x_{i,s}\left(2D_{i}\right)\right]^{T}$
be the local value history vector of $v_{i}$ at $s$ instance of
DAC. First, concatenating $y_{i,s}$ and $y_{j,s},v_{j}\in{\cal N}_{i}$
will obtain
\begin{equation}
\mathbf{q}_{i,s}=\left[\mathbf{y}_{i,s}^{T},\mathbf{y}_{j_{1},s}^{T},\mathbf{y}_{j_{2},s}^{T},\ldots,\mathbf{y}_{j_{\left|{\cal N}_{i}\right|},s}^{T}\right]^{T}.
\end{equation}
Second, concatenating $\mathbf{q}_{i,s}$ will obtain
\begin{equation}
\tilde{\mathbf{q}}_{i}=\left[\mathbf{q}_{i,1}^{T},\mathbf{q}_{i,2}^{T},\ldots,\mathbf{q}_{i,N}^{T}\right]^{T}.
\end{equation}
The vector $\tilde{\mathbf{q}}_{i}$ is constructed by this way so
that each  local value is one time step later than the first  local
value in each row of matrix $\tilde{M}_{i}$. 

Therefore, the coefficient vector can be obtained by inverting the
new  matrix $\tilde{M}_{i}$
\begin{equation}
\mathbf{a}_{i}=-\tilde{M}_{i}^{+}\tilde{\mathbf{q}}_{i}\label{eq:expanded Toeplitz Eq.}
\end{equation}
where $^{+}$ denotes the Moore-Penrose pseudoinverse. Simulation
result indicates that,  the more rows in $\tilde{M}_{i}$, the more
accurate the solution could be. However, increasing the height of
$\tilde{M}_{i}$ after a  limit $n\left|{\cal N}_{i}\right|$ can
not obtain a more accurate result. 


\section{\label{sec:Algorithm-Complexity}Analysis of Algorithm Complexity}

This section is to compare the algorithm complexity of the existing
and proposed distributed optimization  for DAC. 

To solve the problem \prettyref{eq: Opt. FO-DAC Problem}, Xiao \cite{Xiao2004}
proposed two centralized methods: interior-point method and subgradient
method to minimize $\rho\left(W-\mathbf{11}^{T}/n\right)$. Let $n$
be the network size and $m$ be the number of edges. The interior-point
method is iterative and usually finds the optimal solution in $20n\sim80n$
step, at a cost of $(10/3)n^{3}+(1/3)m^{3}$ flops per steps. 

Compared to interior-point method, the subgradient method can be computed
locally but relatively slow. \cite{Boyd2006} proposed a distributed
subgradient method to minimize the sub-dominate eigenvalue of a probability
matrix. This method can be used to optimize the first order DAC after
a little modification. The modified method is given in Algorithm \prettyref{alg:Distributed-FO-DAC-Optimization},
where the function $\mbox{Avg}\left(u^{\left(s\right)}\right)=\frac{1}{n}\sum_{i=1}^{n}\left(u_{i}^{\left(s\right)}\right)$
is implemented by DAC. Given the error tolerance $\epsilon_{ave}$,
the DAC algorithm has to be iterated at least $T_{ave}=O\left(n^{2}log\left(\frac{1}{\epsilon_{ave}}\right)\right)$
times \cite{Zhou2009}, so that the local value vector can converge
to the range $\left\Vert \mathbf{x}\left(T_{ave}\right)-\mathbf{\bar{x}}\right\Vert \leq\epsilon_{ave}$.
The second loop to calculate sub-dominant eigenvector $u_{r}$ is
also similar to the first iteration. It needs to be executed for a
number of times to obtain a result within the error tolerance $\epsilon_{sub}$
\cite{Boyd2006}. Furthermore, the third iteration, which is the subgradient
algorithm,  takes enormous steps to converge and there is no simple
stopping criterion to guarantee a certain level of optimality \cite{Xiao2004}.
Therefore, with a conservative estimation, the times of matrix iteration
might be more than $O\left(n^{4}\right)$.

\begin{algorithm}
\begin{itemize}
\item \textbf{Initialization: }Initialize vector $\mathbf{w}$ with some
feasible entries, for example the maximum degree weights. 
\item \textbf{Repeat} for $t\geq1$

\begin{itemize}
\item $W=I-Q\mathbf{w}^{\left(t\right)}Q^{T}$. 
\item \textbf{Repeat }for\textbf{ $s\geq1$ }to find subgradient $g^{\left(t\right)}\in\mathbb{R}^{\left|{\cal E}\right|}$

\begin{itemize}
\item $u^{\left(s+1\right)}=Wu^{\left(s\right)}$
\item $u^{\left(s+1\right)}=u^{\left(s+1\right)}-\mbox{Avg}\left(u_{i}^{\left(s+1\right)}\right)$
\item $u^{\left(s+1\right)}=u^{\left(s+1\right)}/\left\Vert u^{\left(s+1\right)}\right\Vert $,$s=s+1$
\end{itemize}
\item \textbf{Until } \textbf{$\left\Vert u^{\left(s\right)}-u_{r}\right\Vert _{2}\leq\epsilon_{sub}$}. 
\item $g_{l}=\begin{cases}
-\left(u_{i}-u_{j}\right)^{2} & \mbox{if }\rho\left(W-\frac{\mathbf{11}^{T}}{n}\right)=\lambda_{2}\left(W\right)\\
\left(u_{i}-u_{j}\right)^{2} & \mbox{if }\rho\left(W-\frac{\mathbf{11}^{T}}{n}\right)=\lambda_{n}\left(W\right)
\end{cases}$\\
$l\sim\left(i,j\right)\in{\cal E}$
\item $\mathbf{w}^{\left(t+1\right)}=\mathbf{w}^{\left(t\right)}-\beta_{k}\frac{g^{\left(t\right)}}{\left\Vert g^{\left(t\right)}\right\Vert }$,$t=t+1$.
\end{itemize}
\end{itemize}
\caption{\label{alg:Distributed-FO-DAC-Optimization}Distributively find the
Optimal Matrix for FO-DAC.}
\end{algorithm}


Compare to the enormous number of matrix iterations in Algorithm \ref{alg:Distributed-FO-DAC-Optimization},
there is a significant time reduction by the proposed optimization
if problem \prettyref{eq: Opt. FO-DAC Problem} reduces to find the
best constant. When higher accuracy is needed, we execute $n$ instances
of CFO-DAC and the number of matrix iterations is in the order of
$n^{2}$. 



 Furthermore, if the network topology doesn't change after the optimization,
the estimated eigenvalues  can drive to an estimated  solution $\left(\hat{\epsilon},\hat{\gamma}\right)$
for HO-DAC algorithm, which is faster than the optimal FO-DAC. 


\section{\label{sec:Simulation-of-Eigenvalue}Simulation of Eigenvalue Estimation}

The simulation is taken by the following steps. First, $n$ nodes
are uniformly distributed in an unit square and a link is established
 between any two nodes if their distance is less than a   threshold
$r$. To ensure the generated graph is connected with high possibility,
 $r$ is chosen as $\sqrt{\log_{10}\left(n\right)/n}$ \cite{Li2010}.
Besides, assume  each link is symmetric so that the network graph
is undirected. 

Second, on the random generated network, one or more instances of
CFO-DAC  are executed. Local values obtained in each DAC instance
are stored in local memory of each node. The step length is a constant
shared by all nodes and chosen in the convergence range.

Third,  once sufficient number of local values are obtained, the eigenvalue
estimation algorithm is executed and local Laplacian spectrum is obtained
at each node.

Finally, the performance of eigenvalue estimation is evaluated by
the estimation errors. Before that, each Laplacian eigenvalue $\lambda_{j}\left(L\right)$
is matched with only one eigenvalue $\hat{\lambda}_{k}^{\left(i\right)}\left(L\right)$,
if the distance $\left|\hat{\lambda}_{k}^{\left(i\right)}\left(L\right)-\lambda_{j}\left(L\right)\right|$
is the minimum for all eigenvalues in the estimated local Laplacian
spectrum. Let $e_{i,j}=\min_{l=1,\ldots,D_{i}}\left|\hat{\lambda}_{l}^{\left(i\right)}\left(L\right)-\lambda_{j}\left(L\right)\right|$
be the minimum distance. The mean estimation error of $\lambda_{j}\left(L\right)$
is $\frac{1}{n}\sum_{i=1,\ldots n}e_{i,j}$. 

In Fig.\ref{fig:Box-plot N10 1DAC}, we use the box plot to graphically
illustrate the performance of eigenvalue estimation. The distribution
of log mean estimation errors are obtained from 1000 simulations.
 

\begin{figure*}
\hfill{}\subfloat[\label{fig:10-nodes,-1}10 nodes, 1 DAC instance, with neighbor local
values]{\includegraphics[width=7cm]{\string"D:/Dropbox/PaperWork/Autonomous Opt. DAC/Graph/N10_DAC1_NeiLV\string".pdf}

}\hfill{}\subfloat[\label{fig:10-nodes,-1-1}10 nodes, 1 DAC instance, no neighbor local
values]{\includegraphics[width=7cm]{\string"D:/Dropbox/PaperWork/Autonomous Opt. DAC/Graph/N10_DAC1_NoNeiLV\string".pdf}

}\hfill{}

\hfill{}\subfloat[\label{fig:10-nodes,-2}10 nodes, 2 DAC instance, with neighbor local
values]{\includegraphics[width=7cm]{\string"D:/Dropbox/PaperWork/Autonomous Opt. DAC/Graph/N10_DAC2_NeiLV\string".pdf}

}\hfill{}\subfloat[\label{fig:18-nodes,-18}18 nodes, 18 DAC instance, with neighbor
local values]{\includegraphics[width=7cm]{\string"D:/Dropbox/PaperWork/Autonomous Opt. DAC/Graph/N18_DAC18_NeiLV_grid\string".pdf}

}\hfill{}

\caption{\label{fig:Box-plot N10 1DAC}Box plot of log mean estimation error
of eigenvalues, with/without local value of neighbors. Note that excluding
local values of neighbors will create more outliers as well as increase
the estimation error.}
\end{figure*}


Simulation results show that taking local values from neighbours has
better performance in lower estimation errors. In addition, the estimation
errors will decrease if more instances of CFO-DAC are taken. On the
other hand, the estimation errors  will increase as the network size
becomes larger. However, the numerical errors of $\lambda_{2}\left(L\right)$
and $\lambda_{n}\left(L\right)$ increase very slowly compared to
other eigenvalues in the spectrum. Even their outliers in Fig.\ref{fig:Box-plot N10 1DAC}
have estimated error lower than $1e-8$.

To see how these estimation errors take effect on the performance
of DAC, we conducted another simulation where estimated parameters
$\left(\hat{\epsilon},\hat{\gamma}\right)$  are used to construct
a suboptimal higher-order weight matrix $\hat{H}$. The spectral radius
of $\left(\hat{H}-J\right)$ is plotted and compared with the optimal
one in Fig.\ref{fig:Mean-of-spectral}. 

\begin{figure}
\hfill{}\includegraphics[bb=0bp 0bp 430bp 338bp,width=8.5cm]{\string"D:/Dropbox/PaperWork/Autonomous Opt. DAC/Graph/spectral_radius_Ord123_N8-32\string".pdf}\hfill{}

\caption{\label{fig:Mean-of-spectral}Mean of spectral radius of CFO-DACand
HO-DAC, using optimal parameters and estimated parameters. }
\end{figure}


As shown in Fig.\ref{fig:Mean-of-spectral}, CFO-DAC or SO-DAC using
estimated parameters has similar spectral radius as the one using
optimal parameters. It seems that the numerical errors of $\lambda_{2}\left(L\right)$
and $\lambda_{n}\left(L\right)$ do not decline their performances
dramatically. For third order DAC algorithm, estimated errors of eigenvalues
don't have disastrous impact on the performance. The spectral radius
goes slightly upper than the optimal one after the network size is
larger than 25. It seems that estimating other eigenvalues with low
accuracy is not a critical problem. However, the parameters $\hat{\epsilon}$
and $\hat{\gamma}$ might not be located in the convergence region
if the network size is larger than 32. The simulation of fourth-order
DAC for even larger network up to 40 nodes, reports several divergent
cases.


\section{\label{sub:Conclusion}Conclusion}

In this chapter, we introduced a distributed method to estimate the
optimal parameters for DAC algorithms. However, numerical errors of
these parameters due to quantization  can decline the algorithm performance.
Especially for DAC algorithms with order larger than second, they
are more sensitive to the errors. To mitigate this effect, we introduce
a numerical technique to find the least mean square solution. After
mitigation, the numerical errors of estimated parameters slightly
declines the performance of first order DAC and second order DAC.
For the third order DAC, estimated parameters by local Laplacian spectrum
is still convergent and the algorithm is faster than second order.
 However, if floating point number in double format is used and the
network size is larger than 32 nodes, the numerical errors will be
too large even after mitigation. These findings indicate that the
proposed method is applicable to optimize higher order DAC algorithms
when network size is small.  Otherwise, we should decrease the order
of DAC or increase the accuracy of the floating point number. The
second-order consensus algorithm could be a compromise as it requires
fewer parameters, converges faster than first order DAC and maintains
convergence reliability. In the future, we are intending to investigate
the effect of link failure and other practical aspects while applying
the proposed method to a distributed system. 


